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Abstract — The  Finite  Integration  Technique  (FIT)  ac- 
cording to  T.  Weiland  [1]  is  an  efficient  and  universal 
method  for  solving  a large  scale  of  problems  in  compu- 
tational electrodynamics.  Up  to  now  the  conventional 
formulation  itself  has  had  an  accuracy  order  of  two  with 
respect  to  the  spatial  discretization.  In  this  paper  an 
innovative  extension  to  fourth  or  even  higher  order  is 
presented.  The  convergence  of  the  presented  scheme  is 
demonstrated  by  a general  dispersion  equation  and  sta- 
bility issues  are  discussed.  An  approach  for  a stable  spa- 
tial interface  connecting  regions  of  higher  order  with  the 
standard  FIT  scheme  is  proposed. 

Keywords — -Finite  Integration  Technique,  FDTD,  higher 
order  modeling,  numerical  dispersion 

I.  Introduction 

During  the  last  years  a lot  of  approaches  towards 
higher  order  spatial  finite  difference  ( FD ) schemes  have 
been  developed.  Improved  spatial  discretization  is  nor- 
mally achieved  by  modification  of  the  discrete  curl  oper- 
ator resulting  in  wideranging  spatial  schemes  for  static 
[2],  transient  [3],  [4]  and  frequency-domain  [5]  problems. 
Approaches  using  interpolating  functions  for  fields  with 
cubic  splines  were  also  developed  [6]. 

In  transient  field  analysis  especially  in  the  finite  dif- 
ference time  domain  scheme  (FDTD)  combined  meth- 
ods for  explicit  higher  order  spatial  resolution  and  time 
integration  are  presented  in  [7]- [10].  These  approaches 
use  substitutions  of  higher  order  time  derivatives  by 
spatial  derivatives  leading  to  higher  order  Leap-Frog 
schemes  or  Runge-Kutta  integration  methods.  Recent 
approaches  mix  higher  order  spatial  and  temporal  dif- 
ferencing schemes  to  obtain  a full  fourth  order  accurate 
scheme  for  transient  field  simulation  [12].  A completely 
different  approach  utilizes  multi-resolution  functions  and 
wavelets  for  representation  of  fields  leading  to  higher  or- 
der formulations  in  space  [13].  This  approach  is  currently 
discussed  and  modified  by  various  authors. 

In  this  paper  an  efficient  spatial  formulation  of  arbi- 
trary order  for  the  Finite  Integration  Technique  is  pre- 
sented and  its  applicability  in  the  case  of  a fourth  order 
scheme  (FIT- 4)  is  demonstrated. 

As  one  of  the  key  points  in  the  theory  of  EJT,  the 
modeling  procedure  of  Maxwell’s  Equations  can  be  sep- 
arated into  two  steps.  In  the  first  step,  the  discretiza- 


tion of  the  equations  themselves,  the  so  called  Maxwell’s 
Grid  Equations  are  derived.  Based  on  the  concept  of  grid 
voltages  and  grid  fluxes  they  represent  an  exact  transfor- 
mation of  the  continuous  relations  to  grid  space,  as  the 
integrals  used  are  only  specialized  to  a finite  set  of  in- 
tegration paths  (along  edges  of  the  grids)  or  integration 
volumes  (cells  of  the  grids),  respectively. 

Thus,  the  approximations  of  the  method  do  not  come 
into  effect  until  the  material  matrices  are  introduced  (the 
second  modeling  step).  For  the  derivation  of  these  dis- 
crete analogs  of  the  continuous  constitutive  relations,  the 
integral  state  variables  (fluxes  and  voltages)  have  to  be 
retransformed  to  actual  field  components,  as  will  be  ex- 
plained in  more  detail  later.  For  the  conventional  FIT 
[1],  the  transformation  of  flux  into  voltage  quantities  has 
typically  a second  order  accuracy. 

From  this  point  of  view,  the  path  to  an  extension  to 
higher  order  schemes  has  naturally  to  be  the  following. 
Rather  than  introducing  higher  order  schemes  for  the 
differential  operators  curl  and  div  (as  in  the  FD  liter- 
ature [2]-[12]),  or  defining  higher  order  basis  functions 
for  single  cells  (as  in  p-adaptive  finite  element  schemes), 
utilizing  an  increased  number  of  degrees  of  freedom  per 
cell,  ’’only”  the  material  matrices  have  to  be  replaced 
by  suitable  higher  order  operators.  As  an  important 
consequence,  as  long  as  some  basic  requirements  for  the 
new  material  matrices  are  met,  all  the  well-known  con- 
sistency and  conservation  properties  of  FIT  [14]  can  be 
preserved. 

II.  Basic  Concepts 
A.  FI -Technique 

The  formulation  of  the  Finite  Integration  Technique 
proposed  by  T.  Weiland  [1],  [15]  provides  a general  spa- 
tial discretization  scheme  usable  for  different  electromag- 
netic applications  of  arbitrary  geometry,  e.g.  static  and 
quasi-static  problems  or  calculations  in  frequency-  and 
time-domain. 

The  geometry  is  discretized  on  a dual-orthogonal  grid 
set  consisting  of  the  primary  grid  G (with  the  edges  A/, 
the  facets  A A and  the  material  distribution)  and  the 
so  called  dual  grid  G (containing  the  dual  edges  A l and 
dual  facets  A A).  In  contrast  to  the  vectors  of  elementary 
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field  values  e,  d,  h,  b,  the  FIT  deals  with  the  integral 
expressions 


e = 

/ E • ds, 

Jai 

(la) 

h = 

[ H • ds, 

JAl 

(lb) 

3 = 

[ D-dA , 
Jaa 

(lc) 

3 = 

[ B-dA 
Jaa 

(Id) 

which  form  the  components  of  the  vectors  e,  d,  h and  b, 
being  indicated  by  bows.  The  components  of  the  vectors 
of  electric  voltage  e and  magnetic  flux  6 are  located  on 
the  primary  grid  G)  and  the  components  of  the  vectors 
of  electric  flux  d,  electric  current  j and  magnetic  voltage 
fi  on  the  dual  grid  G (see  Fig.  1). 


tency  of  the  conventional  FIT  formulation 

SC  = 0,  (3a) 

SC  = 0,  (3b) 

C - CT  (3c) 

which  reflect  the  properties  of  its  analytical  pendant. 

For  connecting  the  voltage  and  flux  quantities,  the 
constitutive  relations 

e = (4a) 

= (4b) 

e=Ms-ij  (4c) 

with  the  discrete  material  matrices  M£-i,  M^-i  and 
Mk-i  are  introduced.  They  are  responsible  for  the  dis- 
cretization errors  of  the  method  and  thus  are  the  key 
point  of  the  derivations  in  the  following  sections. 


Fig.  1.  Location  of  electric  voltage  'e  on  edges  and  magnetic  flux 
b on  facets  of  primary  grid  G (a)  and  magnetic  voltage  h on 
edges  and  electric  flux  d and  electric  current  j on  facets  of 
dual  grid  G (b). 


III.  Higher  Order  Material  Relation 

As  explained  before  the  MGEs  deal  only  with  the 
topological  relation  of  the  involved  electric  and  magnetic 
quantities.  Therefore  the  pure  application  of  (2)  is  exact 
i.e.  no  discretization  process  is  applied.  The  constitu- 
tive relations  (4)  connect  fluxes  through  facets  of  one 
grid  with  voltages  along  edges  of  the  corresponding  dual 
grid  which  intersect  these  facets  normally.  The  calcu- 
lation of  the  coupling  coefficients  includes  the  metric  of 
the  grid  as  well  as  the  material  distribution. 

The  scheme  connecting  fluxes  with  voltages  has  to  take 
into  account  the  Maxwellian  continuity  law  of  the  tan- 
gential field  strength  and  normal  flux  density  at  material 
boundaries 


The  FIT  formulation  results  in  the  so  called  Maxwell 
Grid  Equations  (M GEs ) 


C8  = -Se- 

(2a) 

d <=? 

Ch  = — d + J + J ext, 

(2b) 

Sd  = q, 

(2c) 

s6  = o, 

(2d) 

whereby  the  curl  matrices  (C,  C)  and  the  source  matri- 
ces (S,  S)  represent  a summation  scheme  for  the  closed 
line  integral  around  each  cell  facet  and  closed  surface 
integral  over  each  cell  volume,  thus  providing  the  topo- 
logical relation  needed  by  Maxwell’s  integral  equations 
applied  to  the  grid  set.  The  numerical  character  of  the 
spatial  operators  [14]  are  vital  for  the  underlying  consis- 


Etan{ T — EtaniT  (5&) 

Ja{t)  - Htan{r  t)  - 1 Htan(r+,t),  (5b) 

EnormiT  — j =::  Enorm{r  -f,£),  (5c) 

Qa(T  ) — Enormia  — ? £)  Enorm(v  -f,  t)  (5d) 

with  the  surface  current  Ja  and  the  surface  charge  qA- 
These  laws  ensure  in  the  case  of  surface  charge  and  sur- 
face current  free  regions  the  continuity  of  the  tangential 
field  strength  and  normal  flux  density  of  the  electric  and 
magnetic  field. 

A.  Conventional  FIT  Material  Relation 

In  the  following,  a dual-orthogonal  grid  set  with  the 
general  coordinates  u,  v and  w is  regarded.  For  simplic- 
ity reasons,  the  flux  to  voltage  transformation  is  consid- 
ered only  for  the  magnetic  field,  the  conversion  mech- 
anism for  the  electric  field  is  straightforward.  Conven- 
tional FIT  calculates  the  coupling  coefficients  of  mag- 
netic flux  to  magnetic  voltage  in  a two  step  process. 
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1.  The  flux  density  is  derived  from  the  flux  through  the 
related  facet  AA  of  the  primary  grid.  The  definition  of 
the  magnetic  flux  through  a cell  facet 

b = [ B ■ <M,  (6) 

J AA 

ends  in  conventional  FIT  in  the  approximation 

b = b ■ A.4  + 0(A/4),  (7) 


whereby  the  local  flux  density  b is  located  at  the  center 
of  the  related  facet. 

2.  The  locally  calculated  flux  density  is  converted  to 
voltage  by  integrating  it  along  the  corresponding  cell 
edge,  resulting  in  the  multiplication  of  the  flux  density 
value  with  a quotient  of  edge  length  and  proper  averaged 
material  value  respecting  (5c) 

h = li71  B ' ds  + / Ih1  B • ds 

JAh/2  JAh/2 

^^b^+^b^  + OAl2-3).  (8) 

The  length  of  the  dual  cell  edge  is  given  by  A l — 
A/i/2  + A/2/2,  A/i  and  A l2  associated  to  the  adjacent 
cells  of  the  primary  grid.  A full  third  order  scheme  is 
guaranteed,  if  the  grid  is  equally  spaced  and  the  mate- 
rials homogeneous,  otherwise  the  local  conversion  order 
decreases  down  to  0(A12). 

This  two  step  process  ends  in  the  following  formula, 
describing  the  generalized  material  coefficients  for  the 
transformation  of  magnetic  flux  into  magnetic  voltage 


Fig.  2.  Relation  of  magnetic  flux  b with  magnetic  voltage  h in 
conventional  FIT. 


Fig.  3.  Relation  of  electric  flux  ^ with  electric  voltage  'e'  in 
conventional  FIT. 


/. 

J Al. 


A/„ 


H ■ ds=  _ . 

A l„  i^av.^^n  J AAn 


- [ B ■ dA+0(Al2"'3 

n J AAn 


hn 


M-1  ln 

V n.n  ' 


(9) 


with 


1 

Fav. 


A ln 


(10) 


The  necessary  metric  information,  material  distribution 
and  location  of  the  flux  and  voltage  values  for  the  mag- 
netic field  is  displayed  in  Fig. 2,  for  the  electric  field  in 
Fig.  3.  The  resulting  material  coefficients  represent  cell 
inductances,  cell  capacities  and  cell  resistances,  respec- 
tively. 


B.  Principles  of  Higher  Order  Spatial  Discretization 

The  new  higher  order  material  modeling  is  a three  step 
process  utilizing  piece-wise  defined  polynomials  and  fol- 
lowdng  the  basic  ideas  of  the  conventional  FIT. 


• Approximation  of  integrals  associated  to  surfaces  (flux) 
or  edges  (voltage)  with  a suitable  localized  higher  or- 
der polynomial  function  describing  field  strength  or  flux 
density  (in  conventional  FIT  they  are  assumed  to  be 
constant  along  the  edges). 

• Conversion  of  the  derived  field  strength-  or  flux  density 
values  locally  in  their  equivalent  flux  density-  or  field 
strength  values  respectively. 

• Interpolation  of  these  field  values  by  another  localized 
higher  order  polynomial  function  enables  the  calculation 
of  the  desired  voltages  respectively  fluxes. 

Note  that  in  contrast  to  other  higher  order  approaches 
using  wider anging  spatial  differential  operators,  this  new 
approach  leaves  C and  S untouched,  thus  the  properties 
of  FIT  (3)  [14]  still  hold  for  this  discretization  technique. 

In  the  following  section,  a general  method  for  deriving 
higher  order  conversion  schemes  for  surface  based  (flux) 
in  edge  based  (voltage)  quantities  is  discussed  and  an 
exemplarily  approach  for  fourth  order  modeling  (FIT- 4 
scheme)  is  presented.  Once  again  for  simplicity  reasons 
the  discussion  is  restricted  to  the  conversion  scheme  for 
the  magnetic  field,  the  construction  of  the  scheme  for  the 
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electric  quantities  or  the  conductivity  current  is  straight- 
forward. As  seen  before  for  the  conventional  scheme, 
a local  decrease  of  the  convergence  rate  is  inflicted  by 
non-equidistant  grid  spacing  or  inhomogenious  material 
distribution. 


C.  Higher  Order  Flux  to  Voltage  Conversion 
C.l  Derivation  of  Flux 


The  new  approach  assumes  a piece-wise  defined  higher 
order  magnetic  field  strength  function  h(u,  v)  describing 
the  normal  field  component  on  a surface  consisting  of 
facets  A(n^my  The  surface  integral  of  the  assumed  func- 
tion multiplied  by  a material  weighting  function  /i(u,t>) 
approximates  the  magnetic  flux  through  the  surface  con- 
sidering (5b)  and  assuming  surface  current  free  regions. 

For  example  a localized  biquadratic  formulation  of  the 
normal  field  strength  function  on  the  facets  AAu,(. 

A -Au,(.+lii)>  AA„J(._1  jV  A AW(..+1)  and  A can 
be  written  as 


Fig.  4.  Relation  of  the  magnetic  field  strength  value  in  in- 
direction located  at  the  intersection  point  of  dual  edge  A 

and  primary  facet  with  the  magnetic  flux  values  1> 

related  to  the  facets  of  the  primary  grid. 


hw(u,v)  = ai  -f  a2  • u + a3  * t;  + a4  • u2  -h  ■ u2  (11) 

with  the  five  unknowns  a\ , a 2 , • * • , 05 . Starting  with  the 
five  flux  values  oW(i  jV^W(i+Uj),  bW{._,jV  tw  +1),  b^^^ 
and  assuming  an  arbitrary  material  distribution  fiw(u,v) 
on  the  facets  we  postulate: 


scheme  to  all  facets  lead  to  the  five-banded  diagonal  ma- 
trix 


h - 


(15) 


/ 

Jaav 


hw{u,  v)  * liw(u,  v)  du  dv,  (12) 


'(fc.O 


with  (k,l)  = {{i,j),{i  + - 1 ,j),(i,j  + 1)}. 

Evaluating  (12)  within  an  equidistant  grid  with  homoge- 
neous material  distribution  (^(u,^)  = jiw)  leads  to  the 
approximated  flux  through  the  inner  facet  AAW{. 


' W( 


■ j)  *“  Fri) 


AT  AT  , A lU  A lV  . 

■ AluAlv  A o 4 • ~ ) A 


12 


+a5-^A^)  +0(A16). 


(13) 


containing  information  about  metric  and  material  distri- 
bution. 

Fig. 4 displays  the  magnetic  fluxes  bw , metric  coeffi- 
cients A/u,  A lv  and  material  values  \iw  needed  to  com- 
pute hW(i  j)  located  at  the  bary center  of  the  inner  facet. 
Note  that  the  material  inside  one  primary  cell  can  be  dis- 
tributed arbitrarily,  only  the  material  function  fiw(u,v) 
needs  to  be  integratable  over  the  considered  facets. 

For  the  standard  discretization  scheme  (homogeneous 
material  distribution  within  each  cell  as  in  most  FDTD 
implementations),  as  well  as  for  advanced  modeling  ap- 
proaches like  triangular  material  fillings  inside  the  cells, 
explicit  formulas  for  the  matrix  coefficients  can  be  de- 
rived. 


So  the  presented  biquadratic  approximation  leads  to  a 
locally  fourth  order  scheme  for  the  normal  field  strength 
value  hw{ 0,0)  = a%  at  the  intersection  point  of  dual  edge 
A lw  and  primary  facet  A Aw  which  is  also  the  barycenter 
of  this  facet. 

The  resulting  5x5  linear  system 


f <ll\ 

^ u'(i  + l, 3) 

= ■ 

G3 

^Wd.i  + 1) 

G4 

V / 

/ 

(14) 


C.2  Local  Conversion  of  Field  Strength  to  Flux  Density 

The  scheme  described  above  engenders  a vector  of  nor- 
mal magnetic  field  strength  quantities  localized  at  the 
intersection  point  of  dual  cell  edges  and  normal  facets. 
Obeying  (5c)  for  surface  charge  free  domains,  the  normal 
components  of  the  magnetic  flux  density  is  continuous  in 
inhomogeneous  material  distribution.  So  for  calculating 
the  magnetic  voltage,  the  field  strength  component  needs 
to  be  transformed  locally  into  a flux  density  value  (see 
Fig.  5)  leading  to  the  diagonal  matrix 

b = M^h.  (16) 


with  the  local  surface  integration  matrix  MA^  which 
can  be  inverted  for  each  facet.  Applying  this  conversion 


The  described  material  relation  does  not  involve  metric 
coefficients,  so  it  is  free  of  any  discretization  error  due  to 
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& 


Fig.  5.  Local  transformation  of  tangential  continuous  magnetic 
field  strength  into  normal  continuous  magnetic  flux  density  at 
the  intersection  point  p of  normal  facet  and  dual  edge. 


w 


grid  refinement,  however  the  locally  material  smoothing 
inflicts  an  additional  grid  independent  modeling  error. 


Fig.  6.  Integration  of  normal  continuous  magnetic  flux  density 
function^,  (ui)  along  the  dual  edge  to  obtain  the  magnetic 
voltage  hW(k)- 


C.3  Integration  of  Voltage 

Starting  from  magnetic  flux  values  on  facets,  we  de- 
rived in  the  normal  direction  to  these  facets  a contin- 
uous flux  density  quantity.  In  order  to  integrate  this 
flux  density  along  the  related  edge  to  a voltage  value, 
an  interpolation  function  bw(w)  is  assumed  coinciding  at 
every  intersection  point  of  cell  edge  and  corresponding 
dual  facet  with  the  former  calculated  normal  flux  density 
hW(k)  • fiW{ky  For  instance  a localized  quadratic  approach 


for  the  flux  density  function  in  ^-direction 

bw(w)  = ci  + c2  ■ w + c3  • w2  (17) 

with  the  interpolation  conditions 

^iu(O)  =:  5 (18a) 

bw(AlW{k))  — hw^k+l)  * ^(*+1)1  (18b) 

bwi  ^^(fc-i))  i)  * Mu'(fc-i)  (18c) 


results  in  a 3 x 3 linear  equation  system  for  each  dual 
edge.  Having  determined  ci,  C2  and  c3,  the  magnetic 
voltage  can  be  integrated  by  use  of  the  corresponding 
material  distribution  function  fiw(w)  (see  Fig.  6) 


s 


W(fc) 


Ci  + C2  • VJ  + C3  ‘W2 
fiw(w) 


dw. 


(19) 


Integrating  (19)  within  an  equidistant  grid  with  ho- 
mogeneous material  distribution  = jxw)  results 

in 


_ 1 A/3 

h„m  = — (cj  • A lw  + c3-—f)  + 0(A15J.  (20) 

(J>w  *-Z 

In  the  case  of  arbitrary  grid  spacing  or  inhomogeneous 
material  the  quadratic  approach  results  in  a local  con- 
vergence rate  of  C9(A/4”'5). 


Applied  to  all  cell  edges,  the  approach  leads  to  a 
tridiagonal  matrix  converting  magnetic  flux  density-  into 
voltage  quantities 

E = MA(- y„b.  (21) 

C.4  Higher  Order  Material  Matrices  for  High  Frequency 
Problems 

Algorithms  for  high  frequency  problems,  which  will  be 
described  later  in  more  detail,  require  numerical  schemes 
for  flux  to  voltage  transformations.  In  the  case  of  fourth 
order  approximation  and  interpolation  functions,  de- 
noted in  (11)  and  (17),  15  flux  values  are  required  to 
calculate  a single  voltage  value  (see  Fig.  7).  The  above 
described  three  step  scheme  results  in  modified  material 
matrices 


e (22a) 

' v ' 

£ = MAf/MMPM(AAM)-'  S (22b) 

v " 

Mri 

and  the  local  error  of  this  flux-voltage  transformation 
is  0(A/3"'5).  In  the  case  of  an  equidistant  grid  and  ho- 
mogeneous material  distribution,  the  1 5-ban ded  material 
matrices  (see  Fig.  8)  are  positive  semi  definite.  In  terms 
of  numerical  efficiency  it  is  advisable  to  store  the  sur- 
face derivation  matrices  and 

and  the  line  integration  matrices  M and  M &i/e  sepa- 
rately, which  saves  nearly  50%  of  CPU-time  and  memory. 


D.  Higher  Order  Voltage  to  Flux  Transformation 

The  material  matrices  converting  edge  based  (volt- 
age) into  surface  based  (flux)  fields  can  be  derived  in 
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Fig.  7.  Local  transformation  of  magnetic  flux  b into  voltage 
quantities  h with  the  presented  fourth  order  scheme.  The 
calculation  of  one  voltage  value  involves  15  flux  values 


Fig.  8.  Structure  of  the  15-banded  M£-i  matrix  of  a cavity  dis- 
cretized with  10  x 10  x 5 cells  and  PEC  boundary. 


a straightforward  way.  The  basic  steps  are  shortly  de- 
noted for  the  voltage-flux  transformation  of  an  electric 
voltage  component  in  direction:  The  line  integral  of  a 
higher  order  flux-density  function  dw(w)  multiplied  with 
the  material  weighting  function  £w(w)  approximates  the 
voltages  along  the  edges  of  a grid  line.  Converting  the 
flux  density  dw  into  a field  strength  quantity  ew  at  the 
intersection  point  of  dual  facet  and  corresponding  edge 
enables  the  construction  of  a higher  order  field  strength 
function  ew(u,  v)  interpolating  the  derived  field  strength 
values.  Once  again  the  surface  integration  of  the  field 


strength  function  multiplied  with  the  material  weighting 
function  £w(u,v)  enables  the  calculation  of  the  flux  2W 
through  the  inner  facet.  For  illustration  purposes  Fig.  9 
displays  the  15  relevant  electric  voltages  for  calculating 
one  electric  flux.  Note  that  in  contrast  to  the  conven- 


Fig.  9.  Local  transformation  of  electric  voltage  'e'  into  electric 
flux  values  H with  the  presented  fourth  order  scheme.  The 
calculation  of  one  flux  value  requires  15  voltage  values. 


tional  scheme  the  voltage-flux  matrix  Me  is  the  physical, 
but  not  the  numerical  inverse  of  the  flux-voltage  matrix 
Me-i,  i.e.  in  general  holds  M”_\  ^ Me. 

E.  Boundary  Conditions 

Since  the  enhanced  flux- voltage  and  voltage-flux  trans- 
formation leads  to  wideranging  material  matrices,  the 
treatment  of  values  at  the  boundary  of  the  calculation 
area  is  more  enhanced  as  in  conventional  FIT.  PEC 
and  PMC  boundary  properties  are  regarded  as  symme- 
try conditions  for  the  normal  and  anti-symmetry  condi- 
tions for  the  corresponding  tangential  field  strength  val- 
ues at  the  boundary.  Thus,  the  described  algorithm  can 
be  applied  straightforwardly  for  this  kind  of  boundary 
conditions  by  choosing  suitable  even  or  odd  higher  order 
describing  functions  at  the  boundaries.  The  incorpora- 
tion of  open  boundary  conditions  like  Mur’s  ABC  [17] 
or  the  popular  PML-ABC  [18]  follows  the  conventional 
technique. 

F.  Numerical  Efficiency 

Applying  the  described  fourth  order  technique,  the 
simple  diagonal  material  matrices  of  second  order  accu- 
racy are  replaced  by  15-banded  matrices,  so  the  storage 
of  the  material  matrix  and  the  matrix- vector  multiplica- 
tion is  drastically  affected.  Since  the  obtained  matrices 
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are  symmetric,  just  half  of  its  components  have  to  be 
stored,  so  the  memory  requirement  for  the  fourth  order 
scheme  rises  from  N values  (in  the  second  order  case)  to 
5iV  values.  The  CPU-time  increases  from  N to  11 N flops 
for  one  matrix- vector  multiplication.  Assuming  that  y/N 
points  can  be  saved  by  maintaining  the  order  of  accuracy, 
the  fourth  order  scheme  is  numerically  cheaper  than  the 
conventional  one  if  a high  accuracy  solution  is  required. 
Comparing  the  new  approach  with  conventional  higher 
order  finite  differences  approach  [7]  (FDA)  reveals,  that 
the  memory  consumption  of  both  schemes  is  equal  but 
computationally  a single  application  of  the  new  spatial 
operator  costs  14  Flops  (floating  point  operations)  per 
field  component  in  contrast  to  11  Flops  for  the  conven- 
tional FD- 4 scheme  thus  leading  to  a 27%  computational 
overhead. 

IV.  Frequency  Domain  Formulation 

In  the  case  of  time  harmonic  problems,  the  electric 
curl-curl  eigenvalue  equation  can  be  derived  from  (2a) 
and  (2b) 

M,-iCMriCe=Ae  (23) 

with  A = oj2  . The  extraction  of  the  lowest  eigenfrequency 
of  a cavity  discretized  with  N cells  requires  to  shift  the 
nearly  N zero  eigenvalues  caused  by  static  modes  to 
higher  eigenvalues  by  imposing  a graddiv  operation  to 
the  curl-curl  equation.  The  Helmholtz  eigenvalue  prob- 
lem [16] 

(MrtCMri  C + STDsSM”_\  )e  = Ae  (24) 


can  be  constructed  for  the  three  components  of  one  cell 
node  resulting  in  a 3 x 3 eigenvalue  problem  with  the 
three  eigenvalues  co?  = Ai.  The  three  eigenvalues  are 
u>i  = 0 (static  modes)  and  the  two-dimensional  space  of 
eigenvectors  with  the  eigenvalues  c and  reflecting 
the  two  possible  polarization  modes  of  the  plane  wave. 
This  scheme  leads  to  the  formulation  of  a generalized 
grid  dispersion  relation  with  the  eigenvalues 

=Me-,g3)MM-,g2)  ^2sin(^^i))  + 

+ M£-,g3)Mfl-.f1)il)  (2sin(^))2  + 

+ Me-, g2)  (2sin(^))2  , (27a) 

"3  =M£-1(2,2)MM-’(3,3)  (2sin(^y^))  + 

+ M£-1g1)M/t-,g3)  (2sin(^))2  + 

+ M.-.ga)MJ1-,®1)  (2sin(^))2.  (27b) 

It  is  evident,  that  the  transformation  of  surface  based 
to  edge  based  integral  values  has  to  be  treated  sim- 
ilarly for  both  quantities  (i.e.  Me-i = 

etc.)  to  ensure  a physical  descrip- 
tion of  the  plane  wave  propagation. 

Necessary  for  physical  consistency  is  the  convergence 
of  (27)  to  the  analytic  angular  frequency 


with  the  scaling  matrix  D5  ensures  an  appropriate  shift 
of  the  static  eigenvalues  resulting  in  ” ghost  modes” 
which  can  easily  be  identified  using  FIT' s consistency 
relation  (3). 

In  the  present  approach,  where  Me~i  is  a non-diagonal 
matrix  and  its  numerical  inverse  can  not  be  triv- 

ially computed,  the  modified  formulation  [20] 

(CMri  CMri  + SrD'sS)d  = Ad  (25) 

for  the  electric  flux  with  the  modified  shifting  matrix  D's 
can  be  used. 

A.  Grid  Dispersion  Relation 

Assuming  an  infinite  equidistant  grid  with  the  pri- 
mary and  dual  cell  edges  Au,  Av,  Aw  and  a homo- 
geneous material  distribution  with  the  values  e and  p, 
we  consider  the  propagation  of  plane  waves.  Defining 
spatial  phase  factors  Tu  ~ e~jfctiAu,  Tv  = e~jkvAv  and 
Tw  = e~iku,Aw , the  local  version  of  the  curl-curl  matrix 


“ (&«  + kl  + k^)  (28) 

with  ku  = k cos^sin^,  kv  = k sin  <j)  sin  9 and  kw  = 
k cos  0,  which  can  be  revealed  by  means  of  Taylor  expan- 
sion of  the  trigonometric  expressions  of  (27)  with  limA-^o 
proofing  the  order  of  convergence. 

Following  the  described  scheme,  a generalized  local  ma- 
terial matrix 


M W = eMri « = /xMp-i  (29) 


can  be  defined.  Hereby  conventional  FIT  uses  coeffi- 
cients 


M 

M 


(0  _ 

M ~ AVAV 
(0  _ 

W AUAV 


M 


(i)  _ Aw 

<3*3)  _ AUA 


V 


(30a) 

(30b) 

(30c) 


C(')M«IcWM(,2ia(,)  = Xd(l) 

(A  £ 


, . for  the  material  matrix  elements  resulting  in  an  order  of 

' ' 0( A2)  for  (27).  The  proposed  fourth  order  approach  is 
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Fig.  10.  Structure  of  the  electric  curl-curl  matrix  (23)  with  second  order  approach  (a)  and  fourth  order  material  matrices  (b)  of  a cavity 
discretized  with  10  x 10  x 5 cells  and  PEC  boundary  condition. 


described  by  the  following  matrix  elements: 


(31c) 

and  the  Taylor  expansion  of  (27)  reveals  an  overall  spa- 
tial convergence  of  0{ A4).  Figure  11a)  shows  the  en- 
hanced dispersive  character  of  the  fourth  order  modeling 
in  contrast  to  the  second  order  approach  of  a time  har- 
monic plane  wave  propagating  transversely  through  an 
ideal  equidistant  homogeneous  infinite  grid.  The  direc- 
tion dependent  phase  error  at  different  spatial  sampling 
rates  is  displayed  in  11  b).  The  convergence  rate  of  (27) 
using  the  generalized  material  matrix  elements  (30)  re- 
spectively (31)  and  also  of  the  conventional  FD- 4 scheme 
is  displayed  in  a)  and  the  direction  dependent  relative  er- 
ror of  the  FIT- 4 and  the  FD- 4 scheme  in  percent  in  b) 
of  Fig.  12.  The  FD- 4 scheme  demonstrates  a slightly 
lower  dispersion  error  in  every  direction,  the  maximal 
relative  error  of  the  FIT- 4 is  approx  0.2%  larger,  which 
is  quite  negligible. 


V.  Spatial  Stability 

In  order  to  obtain  late-time  stability  of  time-domain 
methods,  the  condition  for  the  so-called  spatial  stability 
[19]  of  the  curl-curl  matrix  (23)  must  hold:  It  states  the 
need  of  real-valued  positive  eigenvalues  of  the  curl-curl 
matrix,  which  is  ensured  in  FIT  by  the  rewritten  form 
of  (23) 

(M;1/2TCM71/2T)r  • (M-^CMJ1/2)®'  = Ae'  (32) 

with  e'  :=  mF2®.  Necessary  for  this  transformation  are 
positive  semidefinite  material  matrices  Mri  and  M^-i, 
which  in  conventional  FIT  is  assured  by  diagonal  matri- 
ces with  non-negative  entries. 

In  the  case  of  a non-equidistant  grid  or  inhomogeneous 
material  distribution  the  fourth  order  method  results 
in  non-symmetric  material  matrices  which  can  lead  to 
complex- valued  eigenvalues  of  (23)  and  an  unstable  up- 
date algorithm  for  time-domain  simulation. 

In  order  to  reobtain  a stable  formulation,  the  sym- 
metrization  of  the  material  matrices  is  enforced  resulting 
in  a local  increase  of  discretization  error.  The  underly- 
ing idea  of  the  symmetrization  process  is  an  averaging  of 
metric  primitives  and  material  values.  Figure  13  displays 
metric  coefficients  needed  for  the  calculation  of  two  ad- 
joined magnetic  voltages  hw and  hw( 2),  Fig. 14  shows 
metric  values,  fluxes  and  material  distribution  for  the 
calculation  of  the  adjoint  magnetic  fluxes  bw( and 

b *u/(2,l)  • 

VI.  Coupling  of  Fourth  and  Second  Order 
Spatial  Regions 

The  conventional  FIT  formulation  offers  a wide  va- 
riety of  enhanced  spatial  discretization  techniques  con- 


, 1 / . 2 / k-U.Au  \ • 2 ( kv^ 

1 + 6lSm  {—)+Sm  [~2 
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Fig.  11.  The  dispersion  characteristic  for  a wave  propagating  in  diagonally  direction  (6  = 54.7°,  <j>  ==  45°)  is  displayed  in  a).  Fig.  b) 
shows  the  phase  error  of  knum  of  the  second  (712)  and  fourth  (714)  material  modeling  at  the  sampling  rates  712/4  =:  = 2,3, 5 and 

10  in  dependence  of  <f>. 


Fig.  12.  Convergence  of  the  spatial  dispersion  relation  (27)  to  the  exact  solution  in  (28).  Figure  a)  displays  the  relative  error  of  the 
angular  frequency  Dnum  of  a plane  wave  traveling  diagonally  through  the  grid  using  second  and  fourth  order  material  matrices  and 
also  the  FD- 4 operator  in  dependence  of  the  number  of  grid  steps  per  wavelength  (n  = A/ A).  Fig.  b)  displays  the  relative  phase 
error  in  percent  of  the  FIT-  4 scheme  (lower  part)  and  of  the  FD- 4 scheme  (upper  part)  at  a sampling  rate  ofn  = A/A=4asa 
function  of  the  direction  of  wave  propagation  in  the  grid. 


cerning  the  treatment  of  non-orthogonal  grids  [20] , sub- 
grids [21],  dispersive  [22],  gyrotropic  [23]  or  non-linear 
material  modeling  [24]  and  a lot  of  other  specialized  tech- 
niques. 

As  seen  before,  the  fourth  order  spatial  modeling  is 
superior  to  the  conventional  scheme  concerning  disper- 
sion, accuracy  and  convergence.  Incorporating  all  these 
features  in  the  fourth  order  technique  would  cause  an 
enormous  numerical  effort  and  programming  and  lead 
to  an  unacceptable  overhead  in  the  computational  pro- 
cess. To  circumvent  this  problem,  a low  reflective,  stable 


and  easy  to  handle  subdomain  technique  combining  the 
advantages  of  second  and  fourth  order  material  model- 
ing has  been  developed.  The  underlying  principle  of  the 
spatial  stability,  explained  in  section  V,  requires  a sym- 
metrical treatment  of  the  components  involved  in  the 
flux-voltage  conversion  process  at  the  interface  connect- 
ing the  second-  and  the  fourth  order  region,  thus  provid- 
ing a symmetry  relation  for  the  involved  quantities. 

At  the  interface  connecting  the  second  and  fourth  or- 
der domains  ” mixed-order”  functions  describing  the  tan- 
gential and  normal  fields  are  used  to  ensure  a symmetric 
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Fig.  13.  Coupling  coefficients  for  calculation  of  magnetic  voltages 
from  magnetic  flux  densities.  Stability  is  guaranteed,  if  cn 2 = 
C/21- 


Fig.  14.  Coupling  coefficients  for  calculation  of  magnetic  field 
strength  values  from  magnetic  fluxes.  Stability  is  guaranteed, 

if  CA12  — CA21- 


interaction  of  the  involved  quantities  from  both  spatial 
domains.  Figure  15  displays  the  structure  of  the  Me-i 
matrix  of  a cavity  modeled  with  the  hybrid  scheme. 

VII.  Example 

The  presented  method  of  fourth  order  spatial  dis- 
cretization is  applied  to  a simple  three  dimensional  rect- 
angular cavity  with  the  dimensions  lmxlmxO.Sm  and 
PEC  boundary  conditions.  The  analytical  resonance  fre- 
quency of  the  configuration’s  lowest  mode  is  fTM1U)  = 
frEu 0 = 212.13  MHz.  The  size  of  the  cell  edges  is  var- 
ied from  A/3  down  to  A/ 17. 

A.  Spatial  Convergence 

To  study  the  convergence  behavior  of  the  fourth  order 
material  matrices  due  to  grid  refinement,  the  resonance 
frequency  is  calculated  in  the  frequency  domain  using 
(24)  respectively  (25)  with  refinement  of  the  cell  edges. 
In  the  equidistant  case  a convergence  rate  of  0(A2  04)  for 
the  conventional  method  and  - as  expected  - of  (9(A4  55) 
in  the  fourth  order  case  can  be  observed.  The  conver- 
gence rates  of  the  eigenfrequencies  are  displayed  in  Fig. 
16  a).  In  the  non-equidistant  case,  whereby  the  spa- 
tial step  is  reduced  from  one  edge  to  the  next  by  5%, 
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Fig.  15.  Structure  of  a mixed  second-fourth  order  M£_i  matrix  of 
a cavity  discretized  with  10  x 10  x 10  cells  and  PEC  boundary. 
Half  of  the  calculation  domain  is  discretized  with  the  fourth 
order  approach,  the  other  one  with  conventional  second  order 
scheme. 

the  frequency  domain  analysis  displayed  in  Fig.  17  a) 
demonstrates  a convergence  rate  of  0( A1*87)  for  the  sec- 
ond order  and  0(A4,18)  for  the  fourth  order  formulation. 

B.  Time  Domain  Convergence 

Higher  order  time  domain  formulations  are  called 
(jV,  V)-schemes,  N describing  the  order  of  temporal 
integration  and  X the  order  of  the  spatial  operators. 
Full  fourth  order  explicit  time  domain  schemes  ((4,4) 
schemes)  require  a fourth  order  spatial  scheme  as  well  as 
a fourth  order  time  integration  method  [25].  For  time- 
domain  analysis  of  the  convergence  behavior,  we  impose 
a Gaussian-formed  pulse  with  200MHz  center  frequency 
and  80MHz  bandwidth  stimulating  the  TMi10  mode. 
The  stimulating  dipole  is  located  at  the  center  of  the 
computational  domain.  The  time  step  size  is  chosen  as 
At  = 0.75- A tcourant  , and  the  resonance  frequency  is  ex- 
tracted by  means  of  Discrete  Frequency  Transformation 
(DFT). 

The  numerically  computed  overall  convergence  order 
in  the  case  of  an  equidistant  grid  is  0(A2,7)  for  the 
(2,2)  Leap-Frog  scheme.  The  fourth  order  scheme  with  a 
fourth  order  version  of  the  Leap-Frog  update  equations 
((4,4)-scheme)  exhibits  an  unexpected  high  convergence 
rate  of  0{ A7,9).  A hybrid  (4,4-2)  Leap-Frog  formula- 
tion [25],  with  a reduced  computational  effort  per  time 
step  uses  fourth  and  second  order  material  matrices  and 
demonstrates  an  overall  convergence  of  C?(A4,8)  (see  Fig. 
16  b).  The  same  analysis  with  a non-equidistant  grid 
(Fig.  17  b)  demonstrates  a convergence  rate  of  O (A2,28) 
for  the  (2,2)  and  0(A4  85)  for  the  (4,4)  scheme. 
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Fig.  16.  Frequency-domain  (a)  and  time-domain  (b)  convergence  of  the  lowest  eigenmode  of  the  rectangular  cavity.  The  frequency 
domain  analysis  was  performed  by  calculating  the  eigenvalues  of  the  curl-curl  matrix  (25).  The  time  domain  analysis  shows  the 
convergence  rate  of  the  (2,2),  (4,4)  and  a hybrid  (4,4-2)  FITD  scheme. 


Fig.  17.  Frequency  Domain  (a)  and  Time  Domain  (b)  convergence  of  the  lowest  eigenmode  of  the  rectangular  cavity  discretized  by  a 
non-uniform  grid  where  the  edge  lengths  are  defined  by  Aj+i  = 0.95Ai  . 


The  improvement  of  the  time-domain  convergence  rate 
of  the  (2,2)  scheme  in  comparison  to  the  second  order 
frequency-domain  modeling  results  from  a partial  com- 
pensation of  time-integration  and  spatial  discretization 
error  having  different  signs.  A similar  effect  is  assumed 
to  be  responsible  for  the  extreme  convergence  rate  of  the 
(4,4)  scheme. 

VIII.  Conclusion 

In  this  paper  a general  extension  of  the  FIT-algorithm 
towards  higher  order  spatial  resolution  is  proposed. 

FIT  formulations  for  all  systems  of  coupled  differen- 
tial equations  following  these  properties  can  be  extended 
to  the  presented  higher  order  technique,  exemplarily  the 
Maxwellian  system  is  discussed. 

The  new  scheme  is  based  on  modified  material  matri- 


ces for  the  transformation  of  grid  fluxes  into  grid  voltages 
and  vice  versa,  which  is  the  only  modeling  step  in  FIT 
where  approximations  are  introduced.  Within  these  ma- 
trices, higher  order  piece- wise  defined  polynomials  are 
applied  for  the  interpolation  of  the  field  quantities,  tak- 
ing care  of  all  physical  continuity  relations.  A gener- 
alized grid-dispersion  equation  is  derived  and  analyzed 
to  demonstrate  the  convergence  of  the  fourth  order  ap- 
proach. The  stability  of  the  new  scheme  is  ensured  by 
the  symmetrization  of  the  resulting  material  matrices.  A 
coupling  technique  for  the  interface  of  second  order  and 
fourth  order  domains  is  discussed. 

In  comparison  to  existing  finite  difference  techniques, 
the  presented  scheme  demonstrates  the  same  dispersion 
characteristics  and  nearly  the  same  computational  cost, 
but  all  the  consistency  and  conservation  properties  of 
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FIT  are  remained,  which  is  not  always  guaranteed  by 
conventional  FD  schemes. 

An  analysis  of  the  lowest  eigenmode  of  a simple  con- 
figuration in  frequency  and  time  domain  using  a fourth 
order  Leap-Frog  scheme  verifies  the  superior  convergence 
rate  of  the  fourth  order  formulation.  The  new  approach 
represents  a vital  enhancement  for  the  applicability  of 
the  FIT  method  to  electrically  larger  problems  and  can 
also  be  used  in  existing  codes  for  error  estimation  pur- 
poses. 
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